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triple-point temperature, of a liquid-like layer at the interface between the coexisting 
solid and vapour phases. 
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1. Introduction 

Nowadays, most of the theoretical studies of the phase behaviour of a classical fluid 
are formulated in the language of the density-functional theory (DFT) PUI2|- Within 
such a theory, the crystalline solid is viewed to be like an inhomogeneous system 
with a periodically modulated density profile n(x), whose free energy is obtained 
through the optimization of a density functional F[n] which is built upon the structural 
properties of the fluid. In particular, a successful recipe for F[n] (sometimes called the 
Hohenberg-Kohn-Mermin (HKM) free energy) is a local mapping into the free energy 
of a homogeneous fluid with a suitably chosen effective density n(x), which is related in 
a non-local way to the real density n(x). At variance with n(x), the smoothed density 
n(x) is a slowly-varying function of the position. This general method has been named 
the "weighted-density approximation" (WDA) [31 Ej . 

This scheme has proved to be sufficient for purely hard-core systems and soft- 
repulsive ones. However, it badly fails in the presence of attractive interactions, where 
it may happen that the solid is being mapped onto a fluid with a density falling into the 
condensation gap where, actually, no homogeneous phase is present. Nor it can be of 
help the fundamental- measure approach of Rosenfeld jHj, which was recently extended 
from hard spheres jHUH] to spherically repulsive interactions [7j, but not yet to attractive 
fluids. In these cases, the only available method is lowest-order perturbation theory pQ, 
which is tantamount to split the HKM free energy into the sum of the density functional 
for a reference system (usually, a hard-core system) and a remainder, containing the pair- 
distribution function of the inhomogeneous reference system. A sensible approximation 
for the latter would allow to draw accurate coexistence lines for the system under 
consideration. In the past, a scheme of this sort has been successfully applied by a 
number of authors to the prototypical case being represented by the truncated Lennard- 
Jones fluid [HI El CHI EEH EE] • 

In the present paper, we prove that this method is effective, although with less 
quantitative success, also for lattice problems. As a case-study, we shall focus on a two- 
dimensional (2D) lattice-gas system which is known ^3] to possess three phases with 
the features of vapour, liquid, and solid. It is argued in Ref. [T3| that the existence in 
this model of a further liquid phase, besides the gaseous one, is made possible by the 
relatively long range of the interparticle attraction. We point out that the choice of a 
2D (rather than 3D) system is only aimed at simplifying the forthcoming analysis of the 
interface problem. In this respect, a crucial test for our DFT will be the prediction of 
surface melting, which is actually in the possibility of the DFT, as is proved by the 3D 
continuum theory of Ref. fT] . 

While the general DFT framework on a lattice (i.e., minimum principle for the 
generalized grand potential plus Ornstein-Zernike relation) is already known UK] , 
we are not aware of any single example of the performance of the lattice DFT in the 
WDA for a system with a realistic phase diagram. We believe that testing the degree 
of sophistication of the lattice-DFT method in a rich context, such as that provided by 
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a three-phase lattice-gas model, can be interesting on fairly general grounds, e.g. for 
weighing up the superiority, if any, of the DFT over other available statistical methods 
like the transfer- matrix [13 a or the cluster variational methods [THj . 

This paper is arranged as follows: After an outline, in section 2, of the main contents 
of the lattice DFT, we describe our system and method in section 3, and present our 
results for the bulk of the system in section 4. Next, in section 5, we analyse the structure 
of the interface between two bulk phases, including a demonstration of the phenomenon 
of surface melting. Further remarks and a brief summary of the main results are given 
in the Conclusions. 



2. Lattice density-functional theory 

We first review the lattice analogue of the classical DFT, which was first considered 
by Nieswand et al. in 1993 [T3j. Like the parent theory in the continuum, the lattice 
DFT is meant to provide a general framework for discussing the statistical properties of 
particles living on a regular lattice, in the presence of a site-dependent external potential 
or, even, of a self-sustained spatial inhomogeneity (like the one which, in a simple fluid, 
becomes manifest at freezing). If lattice sites are allowed to be occupied by at most one 
particle, a general Hamiltonian for our problem is H + J2i with 

H = Y,v{\i-j\)c l c 3 . (1) 

i<j 

We call Ci — 0,1 the occupation number of site i in the lattice, while v(\i — j\) and e$ 
are the pair-interaction and the external potential, respectively. 
The grand-canonical partition function 

~ = £exp - e^c) exp (-(3H) (2) 

{c} V i ) 

(with (3 and \i representing the inverse temperature and the chemical potential, 
respectively) is a lattice functional, namely a function of all components = fi — e* 
of a lattice field, which we call the external field. Then, it is a simple matter to show 
that the (number-) density field is given by: 

. . i dz on 

rii = (Ci) = -ww- = , (3) 

which is a functional of the external field (a temperature dependence is also implied). 
Q = — /3 -1 lnH is the grand potential, also a functional of {/x;}. At variance with the 
continuum case, the density field is a partial, rather than a functional, derivative of S. 

The Hohenberg-Kohn-Mermin theorem, which holds also for lattice gases, ensures 
that there is a one-to-one correspondence, within the space of /x-representable densities, 
between the density field and the external field. Thanks to this theorem, the Legendre 
transform of fi[/x] with respect to its functional variable is a well-defined object, which 
generalizes the concept of Helmholtz free energy to inhomogeneous situations: 



F[n] = + J> 



n, 



(4) 
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Another expression for F[n] is obtained by considering 

7r[n] = i exp ^ /ij[^] c ij ex P (5) 

(with H evaluated at /ij[ n D; which represents the grand-canonical probability density 
for the given n, i.e., the one calculated for p^ = pi[n\. It follows that: 



(6) 

W 

In practice, the ideal-gas system (H = 0) is the only lattice system for which the 
computation of the HKM free energy can be carried out explicitly, with the result: 

f3F id [n] =J2[n i \nn l + (l-n t )\n(l-n l )} . (7) 

i 

In the general case, F[n] is written as the sum of the ideal term and an excess 
contribution .F exc [n] which is to be approximated some way. 

A further density functional can be defined from F[n], which is a sort of generalized 
grand potential: 

nM = F\p]-Y,w = T,*\p]( H +^*\p]-Y,»i c i) > ( 8 ) 

i {c} \ P i / 

where a different symbol, p, is used for the density field to stress the fact that no relation 
is implied between p and p, which should thus be regarded as independent functional 
variables. Instead, we reserve the symbol n for the density field derived from p. 

Quip] is actually the inhomogeneous Gibbs-Bogoliubov functional of the classical 
mean-field theory. Hence, a minimum principle holds, saying that f^p] attains its 
minimum value for a density profile which is precisely the one determined by the given 
external field, namely n. Moreover, the minimum Q^n] is nothing but the grand 
potential Q for the given p. The necessary condition for the minimum reads: 

SI -° (9) 

n Pi=ni 

The minimum principle for f2 M [p] is, besides the HKM theorem, the basic tenet of the 
DFT, being at the heart of a broad family of approximate theories of freezing Every 
such a theory starts from a prescription for F[n], which is then used to determine an 
approximate grand potential for the system from which the thermodynamic properties 
are deduced. 

We conclude our general presentation of the lattice DFT with the Ornstein-Zernike 
(OZ) relation. Taken 

QPexc f)r^ 

■^1 = "'"Si- and ®M = lfc < 10 > 

to be the one- and two-point direct correlation function (DCF), respectively, the formal 
solution of Eq. © reads: 

1 . . 

TH= 7 j-fj y . (11) 

1 + exp {-{3pi - c\ >[n]j 
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Upon introducing the reduced pair distribution function (PDF), 

fti = (12) 



and a further function 

it can be shown (Hj that the following relation follows from Eq. (jllj) : 

hij = Cij + X! Cik n khkj , (14) 

k 

which is the lattice OZ relation (hij = — 1 is called the total correlation function). 
Only a further relation between g and would allow one to determine both functions. 
The importance of is two-fold: on one side, its knowledge permits to recover, 
through the OZ relation, the PDF profile. On the other side, most of the popular 
DFT approximations use an expression of F[n] in terms of the DCF of the fluid. 

A number of simplifications occur for a homogeneous system (i.e., one with e = 0). 
Owing to translational symmetry, the one-point DCF is a constant, c\ [n] = c\(p), 
whereas the two-point DCF is a function of i — j only, that is c^[n] = 02(1 — j,p)- 
Furthermore, Eq. (|14j) can be Fourier-transformed ^7j to give h q = "}2 x h x exp(—'iq 



x 



in terms of C q as: 



h q = ° q ~ , (15) 
q l-pC q 

p being the constant value of the density. 

In the next section, we describe a DFT aimed at reconstructing the phase diagram 

of a realistic lattice gas, that is one with a phase diagram containing, besides a solid 

phase, also two different fluid phases, liquid and vapour. In order to speed up the 

discussion, we have decided to confine most of the technicalities to a few appendices. In 

particular, Appendix A illustrates the lattice counterpart of two celebrated, yet simple, 

theories of freezing: the Ramakrishnan-Yussouff (RY) theory and the Tarazona's WDA. 

We suggest the reading of Appendix A before proceeding to the next paragraph. 



3. Model and method 



We shall work with the t345 model of Ref . JH] • This is a triangular-lattice model with 
a hard core covering first- and second-neighbor sites and a pair attraction ranging from 
third- to fifth- neighbor sites (see Fig. 1 of Ref. JB])- The strength of the attraction 
reduces upon increasing the distance from a reference site: whence, a triangular solid 
is stable at high density and low temperature, with a maximum density of p max = 0.25. 
Upon comparing the solid and the vapour grand potentials, one can easily predict the 
zero-temperature value of the chemical potential at coexistence to be p c (T = 0) = 3t>3, 
v 3 < being the pair-potential value at contact. To be specific, we use hereafter the 
same v values that were considered in Ref. [IHj, namely V3 = — 1.5V, V4 = —1.2V, and 
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v 5 = —V, with V > 0. In that paper, a combination of transfer-matrix calculations and 
Monte Carlo (MC) simulations distinctly showed the existence of a narrow temperature 
interval where the increase of //, starting from large negative values, drives the system 
through a couple of sharp (first-order) phase transitions, i.e., vapour-liquid and liquid- 
solid, as is also revealed by the //-evolution of the number- density histogram at fixed 
temperature. For later convenience, two other models are introduced: the t3 model, 
which is the same as t345 but with V4 = v 5 = 0, and the t model, where also v 3 = 
and only the hard-core interaction is present. At variance with the t345 case, the MC 
simulation supports the existence of a unique fluid phase in both the t and t3 models. 

The first step in a typical DFT calculation is the determination of an accurate 
DCF for the homogeneous system. In fact, an approximate F[n] is usually built upon 
this function (see Appendix A). The fluid DCF is the solution to the homogeneous 
OZ relation plus a closure. For 3D hard spheres, the most celebrated closure of all is 
the Percus-Yevick approximation (PYA), which allows an exact determination of the 
DCF ^Hj- For a lattice system, the mean spherical approximation (MSA) is easier to 
implement than the PYA, since it leads to a smaller set of unknown quantities. As 
a matter of fact, serious convergence problems are encountered when trying to solve 
numerically either the MSA or the PYA of the t345 fluid. Instead, no such problems 
occur for the MSA of the t or, even, the t3 model (see details in Appendix B), while the 
PYA of the t3 model is still intractable. As a result, we see us forced to treat the t345 
model perturbatively, as we are going to see in a moment (note that our derivation of 
the perturbation formula, Eq. 1)16)1 below, will be different from that of Evans pQ). 

Let us write v(\i — j\) as v (\i — j\) + Av(\i — where v Q describes a reference 
system, say the t model, and Av is a remainder. We shall prove that, using a subscript 
for quantities pertaining to the t model, one has, at the lowest order in f3: 



Let v \ = v + XAv be a linear path between v and v, with < A < 1. Accordingly, 
we define H\ = Hq + XAH. Let TT\[n] be the grand-canonical probability density of H\ 
under the condition that the external field takes precisely that value, which 
produces a density of n. Next, we define: 



{c} V H > 

to be the HKM free energy relative to H\. 

In the same spirit of classical Zwanzig perturbation theory JH], we derive an 
approximate expression for F[n] starting from the exact formula: 



F[n]=F [n]+J2M\i-j\)(ciC j ) 



(16) 



i<j 




(17) 




(18) 




(19) 
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where (...)_. is an average over ir\[n] and 

A = J2^^- Ci and B = 5>Ai[n]c,. (20) 

Considering that p\i[n] is unknown, some assumption must be made in order to obtain 
F[n]. In particular, if vq is a hard-core interaction, the r.h.s. of Eq. (j!9|) reduces, at the 
lowest order in ft, to (AH) Q , yielding eventually Eq. iJTfijl . 

We note that, in Eq. (JTHJ), (<HCj) Q = ninjg 0ti j[n\ contains the exact, yet unknown, 
reduced PDF of the inhomogeneous t model. Hence, the above equation is useless unless 
one finds a careful prescription for go,ij, which could be possibly different for the fluid 
and solid phases. Before discussing this point further, we go back for a moment to the 
reference system. 

After obtaining the DCF of the homogeneous t system, we use Eq. (jA.3|) to calculate 
the fluid excess free energy per particle f exc {p)- This quantity, which is a monotonously 
increasing function of the density, ceases to be defined at p ~ 0.21, beyond which 
no MSA solution is actually found. However, this density is too small for allowing a 
description (within the WDA) of the very dense solid. Hence, the problem arises as to 
what criterion should be used in order to extrapolate f cxc {p) beyond that limit. This 
problem is discussed in Appendix B, where two different solutions are proposed. Here, 
suffice it to say that there exists a method to prolong the definition of f exc (p) insofar as 
needed, with all regularity requirements fully met. 

We have sketched in Appendix B the details of a simple DFT (the RY theory [TH] ) 
for the freezing of the t model. However, in order to have a good description of the 
reference system, we have tried to do better than the RY theory. In fact, the stability of 
the liquid phase is a matter of a delicate balance between energy and entropy; hence, an 
accurate representation of the solid free energy is an obvious necessity in all cases where a 
liquid phase is expected. Leaving aside Rosenfeld's fundamental-measure theory, whose 
extension to lattice fluids is not immediate (see, however, the recent contribution |20j), 
we have applied the lattice counterpart of the WDA in the version implemented by 
Tarazona [3], which gives rather good results for the hard-sphere system. This theory 
is reviewed in the Appendices A (general) and C (t model). Here, we provide just a few 
details on the method. 

The hypothesis underlying any WDA is an approximation of the excess free energy 
of the system as Yji n if e ^ c {^i)i where the weighted density fii is a non-local functional 
of the density field, given implicitly by fij = J2j n j w {i ~ j^i)- In turn, the weight 
function w(i — j, p) is such that both the density and the DCF of the fluid are recovered 
in the homogeneous limit. In the Tarazona's WDA, the further assumption is made 
that the weight function is a second-order polynomial in the density p. We thus have 
a well-defined algorithm to build up the excess free energy and, eventually, the density 
functional that is used to trace the conditions for fluid-solid coexistence. 

Once the free energy of the reference system is given, we are left with the problem 
of incorporating the attraction Av into the density functional of the t system using 
Eq. (JTBj). We shall distinguish between the fluid and the solid, although this way the 
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HKM functional will only approximately be the same for all phases (this will have some 
harmful consequences for the interface structure, see section 5). While we obviously 
use for the fluid the reduced PDF of the homogeneous t system, as far as the solid is 
concerned we shall make the (apparently bad) approximation ^o.iiN — 1 outside the 
core, which is the same assumption of the mean-field approximation (MFA). In fact, 
we agree with Mederos et al. [TU] that the PDF of the low-temperature solid is trivial, 
since all the structure (which in a fluid is accounted for by the reduced PDF) is already 
present at the level of ni itself [21]. This is rather obvious at T = 0, where #y = 1 at 
the typical distances of the perfect solid, while being undefined elsewhere. For small, 
but non-zero temperatures, a quasi-random (ideal-gas-like) distribution of interstitials 
and vacancies would extend the result g^ ~ 1 to all distances outside the core region. 

A more refined approximation for the attractive interaction would be that of 
Ref . . This theory uses the same prescription for the solid and the fluid, based on the 
use of the compressibility sum rule. However, the implementation of this method is also 
very difficult and much more involved than ours. In particular, the two algorithms for 
minimization that are described in Appendix B do both require the numerical evaluation 
of the density derivatives of the DFT functional, which is indeed a very difficult task 
to accomplish if the recipe of Ref. JU] is followed. For the sake of truth, we have also 
attempted to use the method of Ref. [H] • This relies on two approximations: i) the use of 
Eq. for the t345 fluid; and ii) the decomposition of the DCF of the inhomogeneous 
t345 system as the sum of the analogous function for the t system and a remainder 
Ac2(i — j, p), assumed to be zero inside the core and MSA-like outside this region. As 
a matter of fact, we found no stable liquid phase by this method. 

Going back to our theory, we write the difference in grand potential between the 
triangular solid (whose density field can be parametrized by means of two numbers only, 
see below) and the fluid with equal T and p as the minimum of: 

Att(n A , n B ) = Afi (t) (n A , n B ) + —k B T {3pv 3 n 2 A + (12pv 4 + 6/3v 5 )n A n B 

5 

+ (9f3v 3 + 12pv A + 6pv 5 )n 2 B - 2p 2 z n Pv n g (n, p) 

n=3 

- (pE z n f3v n g (n, p) + ^- X! Znfiv J^ 90 ™^ ) (n A + 3n B - 4p) 1 . (21) 

V n=3 1 n=3 d P / J 

Note that, in the above equation: i) n A and n B are the number densities in the sublattices 
A and B of occupied and unoccupied sites, respectively (see Appendix B); ii) Afiw is 
the functional for the t model, defined at Eq. (jC.lj) : iii) z n is the coordination number 
for the n-th shell, that is z 3 = 6, z 4 = 12, and z 5 = 6; and iv) g (n, p) is the value 
taken by the reduced PDF of the t fluid at the n-th-neighbour separation. Apart from a 
different density dependence in AQ, the machinery needed for calculating the weighted 
densities n A and n B (and their derivatives) from the densities n A and n B remains the 
same as for the t model, illustrated in Appendix C. 
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The equations for n A and n B are then (see, for comparison, Eqs. (|C.2[) ): 



n 



l-p 



exp [dip) + Pf exc (n A ) + n A (3r c '(n A )^ + 3n B f3f xc '(n B ) dl ' 



B 



pY,Zn(3v n g (n, p) + — z n (3v r 

n n 

+ + (120U4 + 6f3v 5 )n B } ; 



dn A 
dg (n,p) 



dn, 



dp 



n 



l-p 



B 



P 



ex V { Cl (p)+f3r c (n B ) + -n A pf e 



dn, 



dn 



+ n B pr c '(n B 



dn 



B 



B 



dn 



B 



p 

pY,Zn(3v n g (n, p) + — z n pv n 



dg (n, p) 
dp 



+ (A(3v 4 + 2f3v 5 )n A + (Q(3v 3 + 8{3v 4 + Af3v b )n B } . (22) 

In Appendix B, we have outlined two different numerical algorithms for solving the 
minimum problem for a functional of the kind of (|21|). 

We conclude our survey of the method with a few words about the liquid-vapour 
phase transition in the t345 model. The generalized grand potential of the homogeneous 
t345 system is fi^(p) = F(p) — Npp = N(a(p) — pp), where 

f3a(p) = plnp + (1 - p) ln(l - p) + p/?f xc (p) + \p 2 £ z n (3v n g (n, p) . (23) 

n=3 

At low enough temperature, there exists an interval of p values where the minima 
of VLn(p) are in fact two, corresponding to the competing vapour and liquid phases 
(while the deeper minimum yields the physical solution, the other is associated with a 
metastable state). In particular, if we call p v and p\ the related densities, the coexistence 
of the two phases occurs when the minima are equal: 



%(T, p v ) = ty,(T, pi) and Q! AT, p v ) = Q! (T, pi) = . 



(24) 



The above equations are easily identified with the thermodynamic conditions for phase 
coexistence, i.e., equal values of T,P (the pressure), and p for the two phases. This 
will automatically give rise to the Maxwell construction for the pressure and will also 
provide the right position where to cut the non- monotonous profile a'(p) of the chemical 
potential ClS db function of the density. 



4. DFT results: bulk 



In this section, we present the results that we have obtained for the bulk properties of 
the t345 model by the lattice-DFT method outlined in the previous section. In order to 
check them, we have resorted to the MC simulation. In a typical grand-canonical MC 
experiment, a lattice-gas system is driven to equilibrium by a series of moves (creation 
or annihilation of one particle at a time), which are designed in such a way as to 
satisfy detailed balance (for more details, the reader is referred to ^H])- In particular, 
a first-order transition is located at those values of T and p where the number-density 
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hystogram of a large system sample shows two peaks of equal height, signalling that two 
distinct phases are equally stable. 

We first review our results for the t model. We have formulated two different DFTs 
for the freezing of this model, i.e., the RY theory and the WDA of Tarazona. While the 
results of the former are discussed in Appendix B, an outline of the latter can be found 
in Appendix C. Both theories rely on a MSA description for the fluid. Within the WDA 
theory, the densities of the coexisting fluid and solid are found to be pf = 0.1335 and 
p s = 0.1686. Whence, a considerably larger density jump is predicted at the transition 
than given by the RY theory. Anyway, these numbers are still very far from those 
obtained by MC, i.e., pf = 0.172(1) and p s = 0.188(1), indicating that the instability 
of the t fluid against the solid is strongly anticipated in the WDA. As for the chemical- 
potential value at coexistence, the agreement with MC is also poor: while the WDA 
gives (through Eqs. |A~4jl and (fA~5jt ) p c = 1.2655 V, MC yields instead p c = 1.725(5) V. 

In Fig. 1, the local and the weighted density of the t model are separately plotted for 
the two sublattices as a function of the chemical potential. In particular, the weighted 
density takes its larger value in the interstitial region, that is in the B sublattice. This is 
a counterintuitive effect which, however, is not peculiar to the lattice, being also found 
in the continuum (see, for instance, Ref. [I]). 

Moving to the t345 model, we have first checked the existence of two distinct fluid 
phases at low temperature. The liquid-vapour coexistence line is drawn by solving 
Eqs. (J2HJ) and (}2"2|) (see the following Fig. 2). This gives a critical point at t CT = 1.27(1) 
and p cr = 0.079(2) (hereafter, reduced units t = kT/V are used for the temperature). 
Also shown in Fig. 2 is the coexistence line as predicted by the MFA. The latter also 
uses Eq. (J23|) . but with a 1 in place of go(n,p). 

Finally, we have minimized the density functional (J2"T|) in order to obtain the freezing 
and melting lines of the t345 model. It is right at this point that the choice between El 
and E2 (for extrapolating / exc (p) beyond p = 0.21, see Appendix B) becomes crucial. 
In fact, while the solid phase never becomes stable - below a certain temperature - if 
El is adopted, we never run into troubles if extrapolation E2 is used. Anyway, E2 gives 
practically the same results as El at high temperature. 

The complete DFT phase diagram of the t345 model is plotted in Fig. 2 (open 
circles), together with the results of the MFA (crosses) and MC simulation (asterisks). 
To our delight, a triple point eventually shows up in the t345 phase diagram, at 
t tr = 1.145(5) and p tv = 0.122(1), as long as different forms of the perturbation part are 
used in the HKM functional for the solid and for the fluid. In other words, the use of go 
for the description of the reduced PDF of the fluid turns out to be essential for obtaining 
a liquid region in the phase diagram. The liquid phase is unstable if the MFA is used 
also for the fluid. However, the agreement of our DFT with the MC results is mainly 
qualitative: the exact coordinates of the triple point are very different, t^ c = 0.87(1) 
and p^ c = 0.191(1); only the ratio of t tr to t cr is similar. 

As a matter of example, we have plotted in Fig. 3 the /i-evolution of the DFT 
number density at t = 1.2, upon going across the two phase transitions. Finally, Fig. 4 
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shows the DFT phase diagram of the t345 model in the T-\i plane, where we recognize 
the typical fork with two teeth of different length. In the same picture, the MC data 
points of Fig. 2 are also reported for comparison. 

We have studied the t3 model with the same DFT described above in order to check 
the internal consistency of our method. For this model, Eq. ()23j) with v 4 = v 5 = does 
never produce two distinct fluid phases, and the freezing and melting lines are similar 
to those found by the simpler RY theory (see Fig. 5). This result can be rationalized as 
follows: in the t345 model, the existence of attractive sites at the "interstitial" distances 
r 4 and r 5 causes the upper stability threshold of the fluid phase to move up in density 
with respect to the t3 model, thus contributing to unveil the triple point. This effect is 
missing in the t3 model, which thus fails to become a liquid. The conclusion, in perfect 
agreement with MC, is that no liquid phase is present in the t3 model. 

Finally, we make a comment on the possible causes of the quantitative failure of our 
DFT for the t345 model. On one side, one generally expects that mean-field theory works 
well in 3D, less in 2D. One should also not forget that the perturbation formula (fTT)j) is a 
high-temperature approximation and that, at variance with the continuum case, there is 
no Barker-Henderson criterion which can be called for optimizing the hard-core diameter 
of the reference system. On the other hand, also the low quality of the MSA for the 
reference t fluid is partly responsible for the wrong position of the freezing and melting 
lines. To overcome this problem, we have made an attempt of replacing the MSA with 
the hypernetted-chain approximation (HNCA) as a closure for the OZ relation of the 
homogeneous t system. For this model, the HNCA assumes: i) h(0) = h(l) = h(2) = — 1 
(here, the argument is the shell number); ii) C(i — j, p) = h(i — j, p) — In [1 + h(i — j, p)], 
outside the core. In practice, we should also assume that C and h are exactly zero beyond 
a certain distance, and we have chosen this to be the distance of the 38th neighbors 
(i.e., 6\/3). The solution method is iterative: at a given p, we make an initial estimate 
of (7(0), C(l), C(2) and h(3), h(4), . . ., which are then updated using the inverse of the 
Fourier transform ()15|) . Unfortunately, however, this works only up to p = 0.11, which is 
too small a density for allowing us to build an accurate reference-fluid free energy. Just 
in order to appreciate the difference between the two OZ closures, we have plotted in 
Fig. 6 the profiles, for p = 0.1, of the reduced PDF of the t model as given by the MSA 
and by the HNCA, respectively. The comparison with the "exact" MC profile at the 
same density reveals the superiority of the MSA over the HNCA, which overestimates 
the structure of the PDF. However, a good fluid structure is not necessarily accompanied 
by good thermodynamic properties, and this is actually the case of the MSA of the t 
model. 

5. DFT results: interfaces 

Now that we have an accurate density functional for the bulk of the t345 system, 
we move on to consider the structure of the interface between two coexisting bulk 
phases. Many similar calculations have been carried out in the past (see, for instance, 
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Refs. |22[ ITU ITo] ) and, in fact, the development of more and more careful DFT-based 
microscopic descriptions of the density profile across an interface has been historically 
a recurrent leitmotiv j2j- 

Here, two cases are analysed which will deserve a rather different treatment: 
the liquid- vapour interface, i.e., the interface between two homogeneous phases, and 
the solid-fluid interface, which instead separates a broken-symmetry phase from a 
homogeneous one. 

5.1. Liquid-vapour coexistence 

As a first example, we have studied the interface between the coexisting liquid and 
vapour phases of the t345 model. This interface is assumed to lie perpendicularly to the 
y direction. Horizontal layers are labelled with an integer A, which is taken to be zero 
at the "centre" of the interface. Since both phases are homogeneous, the density will be 
uniform along the x direction, its value being a constant, p\, for all sites i of the A-th 
layer. Let p\ and p v be the densities of the coexisting phases at a given temperature 
T. Then, the common value p of the chemical potential is a\p w ) = a'(Pi) (with a(p) 
defined at Eq. (J22J))- For these T and p, the grand potential per site of the bulk vapour 
or liquid is a(p v ) — pp v = a(pi) — pp\. Given that, the generalized grand potential of the 
inhomogeneous system is: 

0fi>] = N X Y, [px In Pa + (1 - Pa) ln(l - p A ) + PxPf exc (px)} 
x 

+ ^Epa £ n J f3Av(\z-j\)g U-j,^^)-N x (3pY,Px, (25) 
1 x j\ieX v J x 

a functional of {p\} being subject to the conditions p\ — > p\ for A — > — oo and p\ — > p v 
for A — > +oo. As is usual practice [JJ, the go function of the inhomogeneous t system at 
the i — j lattice separation is represented by the fluid PDF as calculated for a density 
which is the arithmetic mean of the local densities in i and j. Finally, the grand-potential 
excess per surface particle due to the interface can be estimated as o~(T) = min n E[n], 
where: 

E[n] = ^{0>]-n^)} • (26) 

The calculation of a, which is nothing but the surface tension of the interface under 
consideration, proceeds in two steps: one first optimizes a simple exp ansatz [JS] and 
then refines the calculation via an unconstrained minimization that is accomplished in 
a way analogous to that followed for the bulk. 

We have chosen, for a demonstration, a temperature of t = 1.15, which is slightly 
above the triple-point temperature. For this case, the shape of the liquid-vapour 
interface is plotted in Fig. 7. In this picture, the dotted curve represents the best exp 
profile, while the continuous line is our final optimization. The surface tension is thus 
found to be o~ = 0.0145. By looking at Fig. 7, it appears that the deviation of the density 
profile from the exponential law is actually minute. 
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5.2. Solid-vapour coexistence 

We have first analysed the structure of the solid-fluid interface in the t3 lattice gas by the 
RY theory, as built over the MSA DCF. To be specific, we consider a linear interface 
running along x. Such an interface breaks the translational symmetry along y, thus 
causing the sublattice densities to vary with y. Only very far from the interface, the 
densities recover the bulk values, being those of the solid, say, far below the interface and 
those of the coexisting fluid far above. The horizontal layers are labelled with an integer 
index, A, which increases upon moving from the solid to the fluid region, being zero at 
the interface. We choose e.g. odd A values for those layers where particles are hosted in 
the T = solid. At variance with the bulk case, we must distinguish three sublattices 
since we generally expect different density values at the interstitial sites pertaining to 
the even and to the odd layers. We call C the sublattice formed by the interstitial 
sites in the odd layers, and B the other. Finally, A is the triangular sublattice which is 
occupied in the T = solid. We note that a C site has two adjacent A sites on the same 
layer. Conversely, the two closest A sites of a B site stay on the (odd) layers which are 
respectively below and above the (even) layer which the B site belongs to. 

In the RY theory, the sublattice densities are drawn from Eq. fill)) with p^ = p (given 
by Eq. (jA.5|) ) and a linear density-functionality is assumed for the one-point DCF: 

c «[ n ] = Cl ( p ) + J2c 2 (i-j,p)( nj - p) , (27) 

j 

where i can belong to A, B, or C. In particular, for the odd values of A we have: 

c (1) (A, A) = ci(p) + c 2 (0, p){n A ,\ - p) + 2c 2 (l, p) (n B ,\-i + n c ,x + n B ,x+i - 3p) 

+ c 2 (2, p) (n c ,A-2 + 2n BjA „i + 2n B ,\+i + n c , x+ 2 - 6p) 

+ 2c 2 (3, p) {n A ,x-2 + n A ,\ + n AtX+2 - 3p) ; 
c (1) (C, A) = cx(p) + c 2 (0, p){n c ,x - p) + 2c 2 (l,p) {n B ,\-i + n AX + n Bt x+i - 3p) 

+ c 2 (2, p) (n AtX -2 + 2n BjA _i + 2n BtX +i + n A)X +2 ~ 6p) 

+ 2c 2 (3, p) (n c ,A-2 + n c ,x + n c ,\ +2 - 3p) , (28) 
while, for the even values of A: 

c (1 \B, A) = ci(p) + c 2 (0, p){n B)X - p) + c 2 (l, p) (n A ,A-i + n c ,\-i + 2n B} x + n AtX +i + n c ,\+i - 6p) 
+ c 2 (2, p) (n B> x-2 + n A ,\-x + n G ,\-i + n AiX +i + n c ,\+i + n B , x+ 2 - 6p) 
+ 2c 2 (3, p) (n Bj A- 2 + n BiX + n B)X+ 2 - 3p) . (29) 

Next, the RY density functional is derived from Eq. (jA.6|l . where it must be noted that 

= E( n A,A-p){cU(A\)- Cl (p)} + N x ^ (n B , x -p){c^(B,X)-c 1 (p)} 

A odd A even 

4l - P) {^(C, A) - Cl (p)} . (30) 

Z A odd 
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As is clear, the final expression of the functional (|A.6|) is rather cumbersome and, 
therefore, we do not specify it here. Hence, we directly move to the numerical results. 

We have considered just one temperature value, t = 1.8036. At this temperature, 
the fluid and solid coexistence densities are pf = 0.1000 and p s = 0.1695, respectively. 
Our slab consisted of 61 layers, from A = —30 to A = 30 (at the boundaries, we have 
set the sublattice densities fixed to the solid values for A < —30 and to the fluid value 
for A > 30). To optimize the interface shape, we proceed in two steps: first, we attempt 
a rough optimization by the simple one-parameter ansatz [15] : 

n A ,x = P + 77 77777 ; A odd 

1 + exp(A/t) 

n B ,x = P + 77 77777 ; A even 

1 + exp(A/t) 

n C ,x = P + : — - 77777 • odd) (31) 
1 + exp(A/t) 

The parameter I is chosen in such a way as to make (jA.6|) as low as possible. With 
that profile as a starting point, we run an iterative procedure, similar to that used for 
the bulk, for the unconstrained minimization of fi M [p] — fi M (pf). In the end, we get the 
density profile shown in Fig. 8 (top). At this temperature, the surface tension, given by 
Eq. (JUI, takes the value a = 0.0740(1). 

Next, we move to the t model, as described by the WDA theory outlined in 
Appendix C. From Eq. (|A.5|) . we obtain the following expression for E[n]: 



A odd 



n A>x In h (1 - n A> x) In — h n c ,x In h (1 - n c ,x) In — 

Pv 1 - Pv Pv 1 - Pv 



+ 2 E 



1 nB > X 1 ( ~\ \ 1 1 n B,X 

rz s , A ln h (1 - n B ,x) In 



X even I Pv 1 Pv J A odd 

+ 2 Cl (p v ) £ (n B , x -Pv)+ E [nA,x(3f cxc (n A x) + nc,x(3r c (nc,x)-2 Pv (3r c (p v )} 

X even A odd 

+ 2 E [n B ,x(3f exc (n B ,x) - p v /3/ cxc (Pv)] - (32) 

A even 

For A e.g. even, the only weighted density that matters is n Bj Xj which is defined in terms 
of all densities as: 

n B,x = E n i w {i ~ 3, n B ,\) , (33) 
j 

where i is any particular site on the A-th layer. For odd values of A, one can analogously 
define n A x and n^x- If we adopt the WDA method of Tarazona, then a result similar 
to Eq. (jA.15|) is obtained, giving n Bj x in terms of the auxiliary quantities 

n kB ,x = ^n j w k (i- j) (with k = 0, 1, 2) . (34) 
j 

The explicit expression of (}3~4*j) obviously requires the careful consideration of lattice 
sites j lying progressively farther from the reference site i. As noted in Appendix B, 
a sum like ()34|) should in practice be truncated after a certain value of \i — j\, and we 
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have chosen this to be the distance of the 20th neighbors. Even so, the final formula 
takes too many lines to be specified here, and is therefore omitted. 

The actual minimization of E[n] proceeds in a way analogous to the bulk case, 
described in Appendix C. However, the formulae for the density derivatives of ua,xi ^b,Xi 
and nc,\ are much more involved for the surface than for the bulk case. The outcome 
for the density profile across the interface is shown in Fig. 8 (bottom). Its shape is very 
similar to that of the t3 model, but the surface tension is much larger, our best result 
being a = 0.3182. 

We have finally considered the solid-vapour interface in the t345 model. In 
particular, we are interested in temperature values that are just below the triple-point 
temperature. In such conditions, and as long as surface melting occurs, a thin liquid-like 
film appears at the interface between the solid and the vapour. The functional £[n] for 
the t345 system is the same as for the hard-core model plus the contribution coming 
from the attractive part of the potential: 

>- dg (n,p v )\ 



£[n]=£«[n]- PvE^n<?o< 



n > Pv) + 7^ E Z n(3Vr, 



dp v 



E (n-A,\ + n c ,x - 2p v ) + 2 ^ (n BjX - Pv 



A odd 



A even 



+ E 

^ A odd 



+ E 

A even 



n A,\ E njPAv(\i-j\) +n c ,x E n jl3Av{\i - j\) - 2pl^2z n /3v n g (n,p v ) 

j\i&A,X j\i&C,X n 



n B ,x E n j(3Av{\i 

j\i€B,X 



J\ 



PvE z nPv n g {n, p v ) 



(35) 



The minimization of (f3*K|) is carried out along the same lines as for the reference t 
model, the only difference being in the novel density functionality of E, not in the way 
the weighted density and its derivatives are calculated from the sublattice densities. 

However, we expect a number of oddities to follow from (|33jl because of the different 
functional forms of the solid and fluid free energies. In particular, the minimization 
of (|33j) cannot produce a density profile which, on the A > side of the interface, 
smoothly drops into the vapour one. In fact, contrary to the cases examined before, 
E[n] does not identically vanish when ni takes the constant value p v (even larger is the 
difference, at the triple point, between fi^(pi) and f2^p(p v ), meaning that the obvious 
prerequisite for observing a genuine surface melting is not met). This mismatch can be 
quantified in terms of the difference between p v and the homogeneous solution p^ to 
(Poo) = &u (Pv)- At £ = 1.14, i.e., just below the triple-point temperature, we find 
p v = 0.0356 and p^ = 0.0288 (the difference being smaller at a lower T). 

A way out of this empasse could be that of imposing p^ as boundary value for 
A — > +oo, while maintaining the form ()35j) for E[n]. Obviously, in order to enforce this 
condition, the initial ansatz must be accordingly modified into: 

n A - Poo 



riA,x = Poo + 



1 + exp(A/7) 



(A odd) 
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n B ,\ = Poo + 



n B ~ p, 



(A even) 



1 + exp(A/7) 



nc,x = Poo + 



n B ~ Poo 



. (A odd) 



(36) 



1 + exp(A/7) 



We are perfectly conscious that the solution proposed here just represents a stratagem 
for making Eq. (}35j) suited to describe also the solid surface. A correct description would 
in fact need a unique functional for all phases. 

For future reference, we plot in Fig. 9 (top) the MC outcome for the x-integrated 
densities of the t345 model in a 60 x 128 slab with periodic boundary conditions along 
x and fixed densities at the y boundary. To be precise, the densities are kept fixed at 
the T = solid and vapour values in the eight layers lying on the extreme left and right 
of the picture. The temperature is t = 0.87, i.e., slightly below the exact triple point, 
whereas the chemical potential has been adjusted in order to attain phase coexistence. 
The occurrence of surface melting in the t345 model is demonstrated by the structure 
of the interface in the central part of the picture, which is compatible with that of a 
"modulated" liquid which strongly feels the underlying crystal ordering. 

In Fig. 9 (bottom), we have plotted the density profile across the solid- vapour 
interface at t = 1.14, as calculated through the minimization of S[n]. From a look 
at this figure, we see that there are a few layers, interposed between the solid and the 
vapour, where the values of the sublattice densities are intermediate between those of 
the coexisting solid and vapour and close, on average, to that of the incoming liquid 
(~ 0.121, at t = 1.15). Interestingly, a further evidence (see Fig. 10) goes in support 
of the surface-melting interpretation, namely the existence of a maximum in ub,x near 
A = 0, and of another, less pronounced, in nc\\- These maxima are neither present in 
the initial profile (}36j) nor occur in the interface profiles of the t and t3 models. Anyway, 
the thickness of the molten layer appears to be strongly underestimated by our DFT as 
compared to MC. Moreover, also the comparison with another DFT theory of surface 
melting [TT] does actually lead us to qualify our results as rather poor. 

We are aware that the use of ad hoc boundary conditions in our DFT treatment 
of the solid-vapour interface may cast some doubts on the general significance of the 
results plotted in Figs. 9 (bottom) and 10. Certainly, we are not allowed to draw 
any reasonable estimate of the surface tension from the calculation we have presented, 
which is quantitatively untenable. Notwithstanding the crudeness of our method, we 
nonetheless think that Figs. 9 (bottom) and 10 do indeed catch the genuine behaviour 
of the t345 system. 

6. Conclusions 

In this paper, we have used the lattice DFT method to analyse the phase behaviour 
of a 2D lattice-gas model (named t345) which exhibits a solid, a liquid, and a vapour 
phase. Particles reside on a triangular lattice: occupation of nearest- and next-nearest- 
neighbor sites of a particle is forbidden, while the pair attraction extends from third to 
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fifth neighbors. 

We have built up an accurate solid structure for the purely hard-core model by 
working with the WDA of Tarazona, while the remaining part of the t345 potential has 
been treated as a mean field. This method is expected to provide good results both 
at very low and at very high temperatures, and to offer a not too bad interpolation in 
between. 

As a matter of fact, our theory passes the crucial test of predicting the existence 
of a liquid phase in the t345 model. In particular, the ratio of the triple to the critical 
temperature is found to reproduce the exact value to within 2%. Another successful 
result is the prediction, in agreement with MC simulation, of the disappearance of the 
liquid phase when the range of the attraction is reduced to embrace third neighbors 
only. The main drawback of the theory is in the estimate of the freezing as well as of 
the melting density which, in the worst case, fall short of the exact values by about 35%. 
This inconvenient should be ascribed, besides to the crudeness of the MFA approach 
(also worsened by the low system dimensionality), also to the low quality of the MSA for 
the hard-core fluid (the HNCA is not viable since it does not converge even at moderate 
densities) . 

Having produced a qualitatively sound bulk theory, we have moved to a description 
of the interface structure in the t345 model. The same functional built up for the bulk 
system has been used to describe the coexistence between the solid and the vapour 
phases. Actually, the use of slightly different functional forms for the generalized grand 
potentials of the solid and of the vapour forces us to introduce a spurious boundary 
condition on the vapour side of the interface. If we allow for this artifice, we do in fact 
observe the appearance, just below the triple-point temperature, of a very thin liquid- 
like layer in between the solid and the vapour, which is the sign of the occurrence of 
surface melting in the system. However, it should be admitted that this little evidence 
is not comparable, as for quality, to e.g. that provided by the 3D continuum DFT of 
Ref. HQ- 
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Appendix A. The lattice DFT of freezing — generalities 

In this appendix, we first derive a general expression for the generalized grand potential 
of an inhomogeneous lattice-gas system; this is then used for formulating a lattice DFT 
of freezing. In particular, we show how to adapt the WDA of Tarazona [H] to a lattice 
problem. 

Let us suppose to know the DCF cf^ of a lattice system and the value of its F cxc [n] 
for a given density profile no- Let n\i = no; + AAri;, with An; = n; — no; and < A < 1. 
It then follows from the former of Eqs. (jlU|) that: 

(3F^[n] = PF exc [n ] - £ An, [' dA c 4 (1) [n A ] . (A.l) 

The functional can be similarly obtained, using the second of Eqs. (jlOj) . as an integral 
of c( 2 \ which eventually yields the exact formula: 

PF exc [n] = /3F CX > ] - E c ? } M Arii - £ An;An,- /' dA t dA' cg } [n A ,] , (A.2) 

where the first two terms on the r.h.s. of (JA.2J) identically vanish when choosing uq = 0. 
Actually, an infinite series of terms is hidden behind the last term of Eq. (|A.2|) . each 
containing an order of the DCF as calculated for n . In practice, one could stop this 
infinite regression at the second order by approximating cL [n\>] with q,- [no], and this 
gives the so called Ramakrishnan-Yussouff (RY) theory [T5] . 

Using Eq. (|A.2|) . the excess free energy per particle of a fluid with density p can be 
generally written as: 

/3/ cxc (p) = -- [ P dp' (p - p')c 2 (0, p') , (A.3) 
p Jo 

whereas pf exc {p) is the excess free energy per site. The function c<i is calculated by 
augmenting the OZ relation with a closure, that is a further relation between the total 
and the direct correlation functions. We also quote the expression of C\. 

c 1 (p) = -f3r c (p)-pf3r c '(p), (A.4) 

from which the chemical potential follows through the Eq. 

(3p = \n-^- Cl (p). (A.5) 
1 - p 

In order to study the coexistence between the solid and the fluid, we must 
require equal values of T and p for both phases. Given Eq. (|A.5|) . the departure 
Afl[n] = fi/Jn] — ^ifj,{p) of the generalized grand potential of the solid from that of 
the fluid can be written as: 



0Ml[n] = £ 



Hi In — + (1 — rii) In 



+d(p) J2^i-p)+PF exc [n]-NpPr c (p) , (A.6) 



p l-p 

being N the total number of lattice sites. Every different choice of F exc [n] defines a 
class of (approximate) DFTs. The simplest choice, yet rarely a quantitatively accurate 
one, is the RY theory, which leads to: 



f3AQ[n] = £ 



rii In — + (1 — rii) In 



p l-p 



9^ C 2(* - J, P)(jli - p){rij - p) . (A.7) 
Z i,3 
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The RY theory already represents a considerable improvement over the ordinary MFA, 
which is tantamount to assume c 2 (i — j,p) = for i — j inside the core region, and 
02(1 —j,p) = — (3v(\i — j\) outside the core. At variance with the MFA, the RY theory 
uses a DCF which is adjusted to fit the homogeneous OZ relation as supplemented with a 
proper closure. For instance, in the mean spherical approximation (MSA), one requires 
that g(i — j,p) = inside the core, while still assuming c 2 (z — j,p) = —f3v{\i — j\) 
outside this region. A further possibility would be the Percus-Yevick approximation 
(PYA), which assumes 02(1, p) = g(i,p)[l — exp((3v(\i\))]. 

A more accurate, non-perturbative expression for F exc [n] is obtained by the so 
called weighted-density approximation (WDA) [31 E] , which amounts to approximating 
the exact Eq. (jA.2|) for n Q = as 

F cx >] » -^wdaH = £ mffa) , (A.8) 

i 

where the weighted density rij is implicitly defined by: 

ni = ^2njw(i- j,ni) . (A. 9) 

j 

The weighted density is required to be constant, i.e., fti = p, for a homogeneous system 
of density p (hence, u>(0, p) = J2j w{i—j, p) — 1 for any z); moreover, the weight function 
w must be such that the DCF of the fluid be recovered in the homogeneous limit: 

d2 pexc 
_ j WDA 



dnidnj 



= c 2 (i-j,p). (A.10) 

=p 

With the above requirements, the approximation obtained for F exc [n] is better than any 
truncated DCF expansion 

Using simple calculus, one can translate Eq. (jA.10|) into a differential equation for 
the Fourier transform of w(i, p): 

- ^c 2 (q,p) = 2r c \p)w{q,p) + P r c "(p)w 2 (q,p) + 2 P r c '(p)w(q, p)w'(q, p) . (A.ll) 

Although Eq. (jA.ll|) can be numerically solved for any q and p [4 a , we here adopt the 
simpler recursive method of Tarazona [3], which considers a series expansion of w(q, p) 
in powers of the density. If we stop at the second order, all we need to determine is 

w(i,p) = w (i) + pwi(i) + p 2 w 2 (i) (A. 12) 

from the knowledge of the lower-order terms in the expansions 

Z cxc (p) = fiP + f2P 2 + fsP 3 + ... and c 2 (*,p) = Xo(*) + PXi(i) + P 2 X2^) + ■ ■ ■ (A.13) 
We notice that Eq. (jA.3|) allows us to express fk+i (k = 0,1,2,.. .) in terms of Xk as 

w '«=-(*TO? Xi(,) ' (A - 14) 
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Upon inserting Eqs. (|A.12|) and (|A.13|) into Eq. (|A.11J) and equating term by term, we 
eventually obtain the general formulae: 

w Q (t) = 



wi{q) 
w 2 (q) 



X 1 (g)+W 2 wo(g) + 2/3/ 2 w 2 (g) 

2/3 h (i + w (q)) 

XM + ^hw {q) + AfihwM + 6/3f 3 w 2 {q) + 8/?/ 2 w (gH(g) + 2f3f l w 1 \q) 

2(3f 1 {l + 2w {q)) 

(A.15) 

Given the weight function, the weighted density is explicitly determined from 
Eq. (|A.9|) in terms of the nj as: 

rii = _ , (A.16) 

1 - n u + ^/(l - n Vi ) 2 - 4:n 0i n 2i 

where n ki = Y.j n jWk(i — j) for k — 0, 1, 2. In practice, one uses Eq. (|A.16|) as part of the 
iterative procedure by which the DFT minimum principle is implemented numerically 
(see details in Appendix B). A practical demonstration of the WDA method will be 
given in Appendix C. 
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Appendix B. RY theory for the t3 model 

We hereby describe in detail how to work out a RY DFT for the t3 model. Our first 
task is to determine the DCF of the homogeneous system. We have chosen to solve 
numerically the OZ relation with the MSA closure. Let a x = x and a y = (\/3/2)y 
be the primitive vectors of the triangular lattice (hereafter, we assume a unit lattice 
constant). Then, the reciprocal-lattice vectors are b x = 2tt x and b y = (47ry / 3/3)y. 
Any sum over Born- Von Karman vectors [IT is written, in the infinite-size limit, as an 
integral over the first Brillouin zone: 

1 ^ ... VS f d 2 g y/3 r dq x M/3 dq y 

dq x r dq' y . 2^3^. 
7r 2ix J-tt 2ti 3 
For the t3 fluid, the MSA assumes: i) C(3) = C2(3) = — (3v 3 and C{n) = for all 
n > 3 shells; ii) h(0) = h(l) = h(2) = —1. From ii), three equations are derived for the 
unknown quantities C(0),C(1), and C(2). For instance, the first of these is obtained 
by plugging the OZ relation (|T5j) into the expression h(0) = (1/N)^2 q h q . After a few 
manipulations, we eventually obtain the following set of equations: 



(2vr) 2 i-Tr * i-Tr y 1 - z 1 f 1 (q x , q y ) - z 2 f 2 (q x , q y ) - z 3 f 3 (q x , q y ) 



^l-p)£7(3) = ^ iJ Td fc £d 9 

-6p 2 C(3) = 2:3 r dg dq ACg^gv) 

' ' ' (2tt) 2 7-7T x 7-tt y 1 - z 1 f 1 (q x , q y ) - z 2 f 2 (q x , q y ) ~ z 3 f 3 (q x , q y ) 



-QP 2 C(3) = 7^5 / difc / dg,- r ^ r , (B.2) 

{2ir) z J-tt J-n 1 - z 1 f 1 (q x , q y ) - z 2 f 2 {q x , q y ) - z 3 f 3 (q x , q y ) 

where z x = 2pC(l)/(l -pC(0)), ^ 2 = 2pC(2)/(l -pC(0)), and z 3 = 2pC(3)/(l - pC(0)) 
are auxiliary unknowns. Moreover, 



fi(lx, %) = cos q x + 2 cos (2^) cos % ' 
f2(q x , %) = cos(2q y ) + 2 cos f-gx) cosg^ 



,2 

/s(?x,?») = cos(2g :r ) + 2cosg a; cos(2g y ) . (B.3) 

For a given p, Eqs. (jB.2|) are to be solved recursively: starting from an estimate of z±, z 2 
and z 3 , these quantities are gradually adjusted until the r.h.s. of Eqs. f|B.2jl becomes 
equal to the quantity on the respective l.h.s. with a tolerance smaller than 10 -8 . 

Once the DCF has been determined, we can use e.g. the RY theory to construct 
the generalized grand potential of the t3 model. We call A the triangular sublattice 
that is occupied in the T = crystal, while B includes the rest of the lattice. Then, 
the independent density variables are the two sublattice occupancies, ha and n B , while 
the solid density is p s = (jia + 3tib)/4. The density functional that is to be minimized 
reads: 



= n A In h (1 - n A ) In h 3 

A p 1 — p 



, n B l-n B 
n B m h ( 1 - n B ) In — 

P 1 ~ P 
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- ~ [(c 2 (0) + 6c 2 (3))(n A - p) 2 + 12(c 2 (l) + c 2 (2))(n A - p){n B - p) 

+ 3(c 2 (0) + 4c 2 (l) + 4c 2 (2) + 6c 2 (3))K " P?] , (B.4) 

having omitted to indicate the p dependence of the c 2 values. If, after minimization, AQ 
happens to be negative, then the solid phase is stable, otherwise the fluid will overcome 
the solid in stability. As a result, the locus of the AQ zeroes allows us to draw the 
fluid-solid coexistence line in the T-p plane. 

In order to find out the minimum of AQ, at least two different strategies can be 
pursued whose efficiency turns out, in fact, to be comparable. The first method is 
to lay down, starting from somewhere in the {n A ,nB} space, a fictitious relaxational 
(steepest-descent) dynamics, i.e., 

dAQ 

n A (t + At) =n A (t) - At- — (t) (B.5) 

on A 

and similarly for tib, where At is a conveniently small number. In the long run, the 
sublattice densities eventually stabilize, and this fact will signal that a minimum of 
AQ has been reached (note that there is always the possibility to get stuck in a local 
minimum; so, one should check the nature of the minimum with different At values and 
initial conditions). 

The other method is to solve, by a self-consistent procedure, the non-linear 
equations for the densities, 

n A x = 1 + i^exp [-(c 2 (0) + 6c 2 (3))K - p) - 6(c 2 (l) + c 2 (2))(n B - p)} ; 
P 

n B 1 = l+ exp [-2(c 2 (l) + c 2 {2)){n A - p) - (c 2 (0) + 4c 2 (l) + 4c 2 (2) + 6c 2 (3))(n B - p)\ . 

(B.6) 

In order to reach a better convergence, we have resorted to a mixing scheme: at the 
fc-th step in the iteration, we use the inverse of the r.h.s. of each Eq. (jB.6|) to obtain a 
trial estimate (denoted by a hat) of the densities at the k + 1-th step. Then, we assume 
n A +1 ^ = (1 — q)n A + qn A (and similarly for Ub), where q is a small positive number. 

The RY freezing and melting lines of the t3 model are showed in Fig. 5 as dotted 
lines in the p-T plane. Since there is only one minimum in the fluid generalized grand 
potential, the t3 system shows, according this theory, two phases only - fluid and 
triangular solid, with a density gap becoming narrower and narrower with increasing 
temperatures. In the same figure, we have marked with arrows the densities of the 
coexisting fluid and solid in the t model, namely Pf — 0.1495 and p s = 0.1600 (see 
below). In fact, the t model can be viewed as the infinite-temperature limit of the t3 
model. 

The MSA equations for the t model can be easily adapted from those of the t3 
model. An important thing to notice is that the iterative procedure by which the MSA 
is solved does usually fail to converge beyond a certain density p up which, for the t model, 
is slightly above 0.21. This is a well-known problem in the field of integral equations of 
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classical fluids which, however, is not particularly dangerous in view of the fact that the 
fluid phase loses its stability against the solid well below p np . This notwithstanding, we 
might have the need to extend, as required by the forthcoming WDA of Appendix C, 
the definition of / exc (p) well beyond p up (and even beyond 0.25). To this end, since the 
only obvious constraint to fulfill is regularity, the possible solutions are many. Following 
the proposal advanced in Ref. for hard disks, we could assume, for instance, the 
( met ast able-) fluid pressure to be exactly given, beyond p = 0.21, as: 

(3P l + a'n + tin 2 + c'rf + d'n 4 ._ . 
= 7i ^2 ' ( BJ ) 

p (i-vr 

where n = (27l^/3/3)p = ap is the packing fraction (corresponding to a hard-core 
diameter of 2), while a',b',c', and d! are free parameters. The excess free energy will 
follow from 



t 

which, through the density derivative of Eq. (|A.5J) . is easily proved to be equivalent to 
Eq. (|A.3|) . Upon inserting Eq. (|B.7|) into Eq. (jB.8|) . we eventually arrive at the following 
analytic form: 

[3f exc {p) = ap + bp 2 + — ^— + d ln(l - ap) , (B.9) 

1 — ap 

with other parameters a, b, c, and d. The latter are fixed by requiring a smooth behaviour 
at p = 0.21. 

The problem with the above extrapolation (called El) is that (|B.9|) blows up to 
infinity for p — 1/a ~ 0.276. This could be a serious inconvenient if one needs to 
calculate f exc (p) beyond 1/a. In this case, we resort to a simpler extrapolation (called 
E2), which merely expresses / exc (p) as a fourth-order polynomial beyond p = 0.21. 
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Appendix C. WD A for the t model 

In the present appendix, we show how to build up a WDA theory for the t model. 

Upon inserting (jA.8|) into (|A.6|) . and specializing to the t model, we readily obtain: 

4pAQ(n A ,n B ) n A I - n A \ n B 1 - n B 

— = n A In h (1 - n A ) In h 3 n B m h (1 - n B ) m 

JM p 1 — p p 1 — p 

+ a(p)(n A + 3n B - Ap) + n A (3r c {n A ) + 3n B [3f exc (n B ) - 4p/5/ exc (p) . 

(C.l) 

If we impose the vanishing of the partial derivatives of (jC.l|) . we get the equations for 
n A and n B : 



n A = 1 H exp 

P 

n B = 1 H exp 



,exc/^ •> dn A atwaslfst ^ dn B 



ci(j>) + Pf exc (nA) + n A pr c '(n A )-^ + 3n B pr c \n B/ 

On A On A 

cx{p) + Pr\n B ) + \n A pr c \n A )^ + n B /5/ exc/ (n B ) d " iJ 



3 dn B "' dn B 

(C.2) 

In the above equations, the weighted densities n A and n B are calculated from Eq. (jA.16|) . 
As for their density derivatives, it follows from the original definition ()A.9j) that 

dn A _ _ xfdfioA _ dn 1A _ 2 dn 2A \ . . 

— = (1 - n 1A - 2n 2A n A ) — + n A — + n A — ) , (C.3) 



A \ Uil A Uli>A un A 

and similarly for other derivatives. In Eq. (|C.3|) . n kA = Y,j n jWk(i — j), with i G A and 
k = 0,1, 2. We thus have, for instance, 
dn kA _ 

= w fc (0) + 6w fc (3) + 6w fc (6) + 6w fc (8) + 12w fe (13) + 6w k (l5) + 6w fe (19) + . . . , 

(C.4) 

where we have used the shell number (rather than the distance) as argument for Wk- 
Obviously, in order to make the whole procedure computationally feasible, the sum 
in Eq. (|C.4|) (and any other sum of the same kind) must be truncated at a certain 
distance, and we have chosen to stop summing beyond the distance (equal to 7) of the 
20th neighbors. This is not a problem, however, since the functions rapidly drop to 
zero when increasing the distance from the reference site. 

We remark that a novel feature emerges in the behaviour of w k (i) right when we 
reach the distance of the 20th neighbors, which is not observed at the smaller distances. 
Two different groups of such neighbors are, in fact, to be distinguished: 6 of them are 
symmetry related, as are the other 12. But the value of w k for a site of the first group 
is different from that calculated for a neighbor of the second group (hence we have a 
Wfc(20a) and a w k (20b)). We should wait until the 33th- neighbor shell (at a distance 
of v^T) to observe this feature repeated again. Hence, notwithstanding the potential 
shows radial symmetry, w k is not spherically symmetric and this is the reason why, on 
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a lattice, particular care must always be paid to distinguish translational from spherical 
symmetry, although the exceptions to radial symmetry are, in a sense, rare [21]. Note 
that the reduced PDF behaves similarly to Wk, i.e., g{i—j, p) is not spherically symmetric 
either. 
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Figure Captions 

Fig. 1 : The t model (MSA+WDA). Left panel: total density (p in the fluid phase, n s 
in the solid phase) vs. reduced chemical potential. The sublattice densities, ua and 
ub (dotted lines), are also plotted in the solid region. Right panel, solid region: 
the weighted densities, ua and n B , vs. reduced chemical potential (we have used 
El for extending the definition of f exc (p) beyond p = 0.21). 

Fig. 2 : Phase diagram of the t345 model, using the t model (MSA+WDA) as a 
reference. Two distinct approximations for the perturbation part are compared 
through the respective phase diagrams: Eq. (|21j) (O) vs - the MFA (x). The 
freezing and melting lines are constructed through the use of El at high T and of 
E2 at low T. For comparison, we show as asterisks some MC data points for a 48 x 48 
lattice (MC averages are taken over 5 • 10 5 equilibrium sweeps; the errors affecting 
these points are of the same size as the symbols). The lines connecting the points 
are just a guide for the eye. The arrows pointing downwards mark the densities 
of the coexisting fluid and solid in the t model, as drawn from MSA+WDA. The 
other arrows mark the MC values for the same quantities. 

Fig. 3 : The t345 model (MSA+WDA+perturbation). The picture shows the (3p- 
evolution of the overall density - p for the fluid and n s = (ua + 3ne)/4 for the 
solid - along the isotherm t = 1.2. In the solid region, the real and weighted 
densities are separately plotted for the two sublattices A and B as dotted lines 
(note that % and ha ar e almost indistinguishable). The dotted vertical lines mark 
the points where the phase transitions occur. 

Fig. 4 : The t345 model (MSA+WDA+perturbation). The DFT phase diagram of the 
t345 model (O) as it appears on the T-p plane. The solid-fluid coexistence line is 
of the El type at high T, and of the E2 type at lower T values. For comparison, 
we have also reported as asterisks the MC data for a 48 x 48 lattice. Straight lines 
are drawn through the symbols as a guide for the eye. The inset shows a zoom on 
the triple-point region. 

Fig. 5 : Phase diagram of the t3 model. Two distinct DFTs are contrasted through 
the respective t3 phase diagrams: one theory uses the t model (MSA+WDA) as 
a reference and the attractive interaction as a perturbation (O; the freezing and 
melting lines are constructed through the use of El at high T and of E2 at low 
T); the other theory is MSA+RY (x). MC data for a 48 x 48 lattice and 2 • 10 5 
equilibrium sweeps are shown as asterisks. The straight lines between the points 
are plotted as a guide for the eye. The densities of the coexisting fluid and solid in 
the t model are marked as downward-pointing arrows (long and short arrows are 
for the MSA+WDA and the RY theory, respectively). The other arrows mark the 
MC estimates. 

Fig. 6 : The homogeneous t model. Two distinct closures of the OZ relation are 
compared through the profile of the reduced PDF at p = 0.1: MSA (Q and dashed 
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line) and HNCA (A and dotted line). At the distances r 2 o = 7 of the 20th neighbors 
and r 33 = y/91 of the 33th neighbors, two symbols are reported for each curve (see 
the discussion following Eq. (|C.4J) ). The full dots are the MC data points for a 
48 x 48 sample at /?// = —0.32 (here, the average density is (q) = 0.09995(1) over 
5 ■ 10 5 equilibrium sweeps). Inset: a magnification of the large-distance region. It 
distinctly appears that the MSA PDF is of an overall better quality than the HNCA 
one. 

Fig. 7 : The t345 model (MSA+WDA+perturbation). The figure shows the density 
profile across the interface between the coexisting liquid and vapour phases at 
t = 1.15. The starting point of the functional minimization is an exp profile (dotted 
line); The open dots (which, for clarity, are joined by a continuous line) are the 
final outcome of the optimization. It turns out that the difference between the two 
curves is very small. 

Fig. 8 : Top: the t3 model (MSA+RY), density profile across the solid-fluid interface 
at t — 1.8036 (here, pf = 0.1000 and p s = 0.1695). Bottom: solid- fluid interface 
in the t model (MSA+WDA). In both panels, the optimal exp profile (dotted line) 
is contrasted with the outcome of an unconstrained E[n] minimization (open dots 
and continuous line). 

Fig. 9 : Top: MC density profile across the solid-vapour interface of a t345 lattice 
system at t = 0.87, i.e., just below the triple-point temperature. To reach 
coexistence, the chemical potential is set equal to \i — —4.479 V. For a simulation 
box of 60 x 128, as many as 2 ■ 10 6 equilibrium sweeps were produced. The dotted 
line marks the average density over couples of adjacent layers. Near the centre of 
the picture, the maximum in the interstitial density (which is the bottom of the 
modulation) is the sign of a liquid-like behaviour. Bottom: DFT results for the t345 
model (MSA+WDA+perturbation). The density profile across the solid- vapour 
interface is shown at t = 1.14: optimal exp profile (dotted line) vs. unconstrained 
S[n] minimization (open dots and continuous line). 

Fig. 10 : The t345 model (MSA+WDA+perturbation). The profile of ub,\ and nc,\ 
vs. A for the same solid- vapour interface being represented in Fig. 9 (bottom). The 
unconstrained minimum-S[n] profile (open dots and continuous lines) is compared 
with the best exp ansatz (dotted line). The maxima near A = are a sign of the 
local onset of a liquid-like behaviour. 
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